/*
 * Copyright 1993-2010 NVIDIA Corporation.  All rights reserved.
 *
 * Please refer to the NVIDIA end user license agreement (EULA) associated
 * with this source code for terms and conditions that govern your use of
 * this software. Any use, reproduction, disclosure, or distribution of
 * this software and related documentation outside the terms of the EULA
 * is strictly prohibited.
 *
 */
 
/*
    Volume rendering sample
    This sample loads a 3D volume from disk and displays it using
    ray marching and 3D textures.
    Note - this is intended to be an example of using 3D textures
    in CUDA, not an optimized volume renderer.
    Changes
    sgg 22/3/2010
    - updated to use texture for display instead of glDrawPixels.
    - changed to render from front-to-back rather than back-to-front.
*/

/*  This file includes modifications to the SDK 
 *  Several volume rendering techniques were added 
 *  Authors: Jim Pelton and Vanessa Gertman
 *  CS464 Fall 2011 Final Project
 *  Prof. Alark Joshi
 *  Boise State University
 *
 */

//#ifdef WIN32
//#include <Windows.h>
//#endif

// OpenGL Graphics includes
#include <GL/glew.h>
#if defined (__APPLE__) || defined(MACOSX)
#include <GLUT/glut.h>
#else
#include <GL/freeglut.h>
#endif

#include "volumeRender_param.h"


// CUDA utilities and system includes
#include <cutil_inline.h>    // includes cuda.h and cuda_runtime_api.h
#include <cutil_gl_inline.h> // includes cuda_gl_interop.h// includes cuda_gl_interop.h
#include <rendercheck_gl.h>

// CUDA Includes
#include <vector_types.h>
#include <vector_functions.h>
#include <driver_functions.h>

// Shared Library Test Functions
#include <shrUtils.h>
#include <shrQATest.h>

#include <iostream>
#include <fstream>
#include <map>
#include <string>

struct IsoPair{
    int isoval;
    float isoalpha;
};

#define MAX(test, old) ((test > old) ? test : old)

//in our special case we could optimize Clamp, because value<min is always false since
//min is 0 but value is an unsigned type.
#define Clamp(value, min, max) (value < min) ? min : (value > max) ? max : value

#define MAX_EPSILON_ERROR 5.00f
#define THRESHOLD         0.30f

typedef unsigned int uint;
typedef unsigned char uchar;


// Define the files that are to be save and the reference images for validation
const char *sOriginal[] =
{
    "volume.ppm",
    NULL
};

const char *sReference[] =
{
    "ref_volume.ppm",
    NULL
};

const char *sSDKsample = "CUDA 3D Volume Render";

//default file if none is given on CL opts.
char *volumeFilename = "bonsai128x128x64.raw";
cudaExtent volumeSize = make_cudaExtent(128, 128, 64);
typedef unsigned char VolumeType;
char *transFuncFileName = "transfer.tf";
int tfLines = 4;


dim3 blocksize(16,16);
dim3 gridSize;
struct vrParams kargs;

int method = 1;
uint width = 512, height = 512;
float3 viewRotation;
float3 viewTranslation = make_float3(0.0, 0.0, -4.0f);
float invViewMatrix[12];

bool linearFiltering = true;

int g_singleDimVolumeDataArrayLength;
float *h_depthBuffer;
float *d_depthBuffer;

vector<IsoPair> *g_isoValues;
vector<float4> *g_transValues;
map<string, float> *g_options;


GLuint pbo = 0;     // OpenGL pixel buffer object
GLuint tex = 0;     // OpenGL texture object
struct cudaGraphicsResource *cuda_pbo_resource; // CUDA Graphics Resource (to transfer PBO)

unsigned int timer = 0;

// Auto-Verification Code
const int frameCheckNumber = 2;
int fpsCount = 0;        // FPS count for averaging
int fpsLimit = 1;        // FPS limit for sampling
int g_Index = 0;
unsigned int frameCount = 0;
unsigned int g_TotalErrors = 0;
bool g_Verify = false;
bool g_bQAReadback = false;
bool g_bQAGLVerify = false;
bool g_bFBODisplay = false;

// CheckFBO/BackBuffer class objects
CheckRender *g_CheckRender = NULL;

int *pArgc;
char **pArgv;

int parseTrans(char *fname, map<string,float> *options, vector<float4> *myTfunc, vector<IsoPair> *myIsoVals);

extern "C" void setTextureFilterMode(bool bLinearFilter);
extern "C" void initCuda(void *h_volume, float *h_grads, float4 *h_transFunc, float *d_dbuf, size_t tfsize, cudaExtent volumeSize);
extern "C" void freeCudaBuffers();
extern "C" void render_kernel(dim3 gridSize, dim3 blockSize, uint *d_output, float *d_dbuf, int method);

extern "C" void copyInvViewMatrix(float *invViewMatrix, size_t sizeofMatrix);
extern "C" void copyKParamsToDevice(struct vrParams *par, size_t sizeofParams);

extern "C" void copyDepthBufferToHost(float *h_dest, float *d_src, size_t sz);


void initPixelBuffer();
void computeGradient(char *, float*);

void AutoQATest()
{
    if (g_CheckRender && g_CheckRender->IsQAReadback()) {
        char temp[256];
        sprintf(temp, "AutoTest: Volume Render");
        glutSetWindowTitle(temp);
	shrQAFinishExit2(false, *pArgc, (const char **)pArgv, QA_PASSED);
    }
}

unsigned int fpsTimer = 2;
float thisTime;
float lastTime;
void computeFPS()
{
	
	frameCount++;
    fpsCount++;
    if (fpsCount == fpsLimit-1) {
        g_Verify = true;
    }
	
    if (fpsCount == fpsLimit) {
        char fps[256];
        float ifps = 1.f / (cutGetAverageTimerValue(timer) / 1000.f);
        sprintf(fps, "%sVolume Render: %3.1f fps", 
                ((g_CheckRender && g_CheckRender->IsQAReadback()) ? "AutoTest: " : ""), ifps);  
        

        
        
        glutSetWindowTitle(fps);
        fpsCount = 0; 
        if (g_CheckRender && !g_CheckRender->IsQAReadback()) fpsLimit = (int)MAX(ifps, 1.f);

        cutilCheckError(cutResetTimer(timer));  
		
		thisTime = cutGetTimerValue(fpsTimer);
		float dt = thisTime-lastTime;
		if (dt >= 1000.0f){
			//fprintf(stdout, "fps: %.2f\n", ifps);
			shrLog("fps: %.2f", ifps);
		}

        AutoQATest();
    }
	
	lastTime = thisTime;
    
}

// render image using CUDA
void render()
{
	copyInvViewMatrix(invViewMatrix, sizeof(float4)*3);


    // map PBO to get CUDA device pointer
    uint *d_output;
	// map PBO to get CUDA device pointer
	cutilSafeCall(cudaGraphicsMapResources(1, &cuda_pbo_resource, 0));
    size_t num_bytes; 
    cutilSafeCall(cudaGraphicsResourceGetMappedPointer((void **)&(d_output), &num_bytes,
						       cuda_pbo_resource));
    //printf("CUDA mapped PBO: May access %ld bytes\n", num_bytes);

    // clear image
    cutilSafeCall(cudaMemset(d_output, 0, width*height*4));

    // call CUDA kernel, writing results to PBO
    render_kernel(gridSize, blocksize, d_output, d_depthBuffer, method);

    cutilCheckMsg("kernel failed");

    cutilSafeCall(cudaGraphicsUnmapResources(1, &cuda_pbo_resource, 0));
}

// display results using OpenGL (called by GLUT)
void display()
{
    cutilCheckError(cutStartTimer(timer));  

    // use OpenGL to build view matrix
    GLfloat modelView[16];
    glMatrixMode(GL_MODELVIEW);
    glPushMatrix();
        glLoadIdentity();
        glRotatef(-viewRotation.x, 1.0, 0.0, 0.0);
        glRotatef(-viewRotation.y, 0.0, 1.0, 0.0);
        glTranslatef(-viewTranslation.x, -viewTranslation.y, -viewTranslation.z);
    glGetFloatv(GL_MODELVIEW_MATRIX, modelView);
    glPopMatrix();

    invViewMatrix[0] = modelView[0]; invViewMatrix[1] = modelView[4]; invViewMatrix[2] = modelView[8]; invViewMatrix[3] = modelView[12];
    invViewMatrix[4] = modelView[1]; invViewMatrix[5] = modelView[5]; invViewMatrix[6] = modelView[9]; invViewMatrix[7] = modelView[13];
    invViewMatrix[8] = modelView[2]; invViewMatrix[9] = modelView[6]; invViewMatrix[10] = modelView[10]; invViewMatrix[11] = modelView[14];

    render();


    // display results
    //glClear(GL_COLOR_BUFFER_BIT);
    glClearColor(1,1,1,1);

    // draw image from PBO
    glDisable(GL_DEPTH_TEST);

    glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
#if 0
    // draw using glDrawPixels (slower)
    glRasterPos2i(0, 0);
    glBindBufferARB(GL_PIXEL_UNPACK_BUFFER_ARB, pbo);
    glDrawPixels(width, height, GL_RGBA, GL_UNSIGNED_BYTE, 0);
    glBindBufferARB(GL_PIXEL_UNPACK_BUFFER_ARB, 0);
#else
    // draw using texture

    // copy from pbo to texture
    glBindBufferARB(GL_PIXEL_UNPACK_BUFFER_ARB, pbo);
    glBindTexture(GL_TEXTURE_2D, tex);
    glTexSubImage2D(GL_TEXTURE_2D, 0, 0, 0, width, height, GL_RGBA, GL_UNSIGNED_BYTE, 0);
    glBindBufferARB(GL_PIXEL_UNPACK_BUFFER_ARB, 0);

    // draw textured quad
    glEnable(GL_TEXTURE_2D);
    glBegin(GL_QUADS);
    glTexCoord2f(0, 0); glVertex2f(0, 0);
    glTexCoord2f(1, 0); glVertex2f(1, 0);
    glTexCoord2f(1, 1); glVertex2f(1, 1);
    glTexCoord2f(0, 1); glVertex2f(0, 1);
    glEnd();

    glDisable(GL_TEXTURE_2D);
    glBindTexture(GL_TEXTURE_2D, 0);
#endif

    if (g_CheckRender && g_CheckRender->IsQAReadback() && g_Verify) {
        // readback for QA testing
        shrLog("\n> (Frame %d) Readback BackBuffer\n", frameCount);
        g_CheckRender->readback( width, height );
        g_CheckRender->savePPM(sOriginal[g_Index], true, NULL);
        if (!g_CheckRender->PPMvsPPM(sOriginal[g_Index], sReference[g_Index], MAX_EPSILON_ERROR, THRESHOLD)) {
            g_TotalErrors++;
        }
        g_Verify = false;
    }

    glutSwapBuffers();
    glutReportErrors();

    cutilCheckError(cutStopTimer(timer));  

    computeFPS();
}

void idle()
{
    glutPostRedisplay();
}

void printStuff()
{
    char *thing =
    	"|M: m=%d                                  |\n"\
    	"|den=%.2f brt=%.2f  trOff=%.2f trSc=%.2f  | \n" \
        "|Bound:   kgc=%.2f  kgs=%.2f   kge=%.2f   |   \n" \
        "|Depth:   kds=%.2f  kde=%.2f              | \n" \
        "|Tone:    ktd=%.2f  ktb=%.2f   kty=%.2f   |   \n" \
        "|MIDA:    mg=%.2f                         | \n" \
        "|Sil:     ksc=%.2f  kss=%.2f   kse=%.2f   | \n" \
        "|Light:   lx=%.2f   ly=%.2f    lz=%.2f    |   \n" \
        "|Material:amb=%.2f  dif=%.2f   spc=%.2f   |   \n" \
        "|Classify:iv1=%.2f  ia1=%.2f              | \n" \
        "|         iv2=%.2f  ia2=%.2f              | \n" \
        "|         iv3=%.2f  ia3=%.2f   fo=%.2f    |   \n" \
        "|\n";
    shrLog(thing, 
        method, kargs.density, kargs.brightness, kargs.transferOffset, kargs.transferScale,
        kargs.kgc, kargs.kgs, kargs.kge, 
        kargs.kds, kargs.kde,
        kargs.ktd, kargs.ktb, kargs.kty,
        kargs.mida_gamma,
        kargs.ksc, kargs.kss, kargs.kse,
        kargs.lightx, kargs.lighty, kargs.lightz, kargs.mat_a, kargs.mat_d, kargs.mat_s,
        kargs.isoValues[0], kargs.isoAlphas[0], kargs.isoValues[1], kargs.isoAlphas[1], 
        kargs.isoValues[2], kargs.isoAlphas[2], kargs.fallOff);
}


void writeDBufToFile()
{
    std::ofstream file("dbuf.txt");  
    if (!file.is_open()){ std::cout << "file wasn't opened" << std::endl; return; }
    int i = 0;
    while (i < volumeSize.width*volumeSize.height*volumeSize.depth)
    {
        file << h_depthBuffer[i++] << std::endl;
    }
    file.close();
}


void keyboard(unsigned char key, int x, int y)
{
    switch(key) {
        case 27:
	    shrQAFinishExit2(false, *pArgc, (const char **)pArgv, QA_PASSED);
            break;
        case '\\':
            linearFiltering = !linearFiltering;
            setTextureFilterMode(linearFiltering);
            break;

        /****\\**************************/
        /* DENSITY/BRIGHTNESS         */
        /******************************/
        case '=':
            kargs.density += 0.01f;
            break;
        case '-':
            kargs.density -= 0.01f;
            break;

        case '_':
            kargs.brightness += 0.1f;
            break;
        case '+':
            kargs.brightness -= 0.1f;
            break;

        /****\\***************************/
        /* TRANSFER OFFSET/SCALE       */
        /*******************************/
        case ']':
            kargs.transferOffset += 0.01f;
            break;
        case '[':
            kargs.transferOffset -= 0.01f;
            break;
            
        case '{':
            kargs.transferScale += 0.01f;
            break;
        case '}':
            kargs.transferScale -= 0.01f;
            break;

        /***\\********************/
        /* EDGE ENHANCEMENT    */ 
        /***********************/
        case 'a':
            kargs.kgc += 0.1f;
            break;
        case 'z':
            kargs.kgc -= 0.1f;
            break;

        case 's':
            kargs.kgs += 0.1f;
            break;
        case 'x':
            kargs.kgs -= 0.1f;
            break;

        case 'd':
            kargs.kge += 0.1f;
            break;
        case 'c':
            kargs.kge -= 0.1f;
            break;
        
        /*****\\********************/
        /* DEPTH COLOR BLENDING  */
        /************************/
        case 'f' :
            kargs.kds += 0.1f;
            break;
        case 'v' :
            kargs.kds -= 0.1f;
            break;

        case 'g' :
            kargs.kde += 0.1f;
            break;
        case 'b' :
            kargs.kde -= 0.1f;
            break;
        
        /**\\************************/
        /* TONE SHADING           */
        /**************************/
        case 'h':
            kargs.ktd += 0.1f;
            break;
        case 'n':
            kargs.ktd -= 0.1f;
            break;

        case 'j':
            kargs.ktb += 0.1f;
            break;
        case 'm':
            kargs.ktb -= 0.1f;
            break;

        case 'k':
            kargs.kty += 0.1f;
            break;
        case ',':
            kargs.kty -= 0.1f;
            break;

        /************************/
        /* MIDA GAMMA           */
        /************************/
        case 'l':
            kargs.mida_gamma += 0.1f;
            break;
        case '.':
            kargs.mida_gamma -= 0.1f;
            break;

         /***\\**********************/
         /* LIGHTING              */
         /*************************/

        case 'q':
          kargs.lightx += 0.1f;
          break;
        case 'w':
          kargs.lightx -= 0.1f;
          break;
        case 'e':
          kargs.lighty += 0.1f;
          break;
        case 'r':
          kargs.lighty -= 0.1f;
            break;
        case 't':
            kargs.lightz += 0.1f;
            break;
        case 'y':
            kargs.lightz -= 0.1f;
            break;

        case 'Q':
            kargs.mat_a += 0.05f;
            break;
        case 'W':
            kargs.mat_a -= 0.05f;
            break;

        case 'E':
            kargs.mat_d += 0.05f;
            break;
        case 'R':
            kargs.mat_d -= 0.05f;
            break;

        case 'T':
            kargs.mat_s += 0.05f;
            break;
        case 'Y':
            kargs.mat_s -= 0.05f;
            break;

        /****\\*********************/
        /* CLASSIFICATION        */
        /*************************/

         case 'o':
         	kargs.fallOff += 0.1f;
         	break;
         case 'p':
         	kargs.fallOff -= 0.1f;
         	break;

         case 'A':
             kargs.isoValues[0] += 1.f;
             break;
         case 'Z':
             kargs.isoValues[0] -= 1.f;
             break;
         case 'S':
             kargs.isoAlphas[0] += 0.05f;
             break;
         case 'X':
             kargs.isoAlphas[0] -= 0.05f;
             break;

         case 'D':
             kargs.isoValues[1] += 1.f;
             break;
         case 'C':
             kargs.isoValues[1] -= 1.f;
             break;
         case 'F':
             kargs.isoAlphas[1] += 0.05f;
             break;
         case 'V':
             kargs.isoAlphas[1] -= 0.05f;
             break;

         case 'G':
             kargs.isoValues[2] += 1.f;
             break;
         case 'B':
             kargs.isoValues[2] -= 1.f;
             break;
         case 'H':
             kargs.isoAlphas[2] += 0.05f;
             break;
         case 'N':
             kargs.isoAlphas[2] -= 0.05f;
             break;
                
        /****\\*******************/
        /* silhouettes         */
        /***********************/
         case 'J':
             kargs.ksc += 0.1f;
             break;
         case 'M':
             kargs.ksc -= 0.1f;
             break;

         case 'K':
             kargs.kss += 0.1f;
             break;
         case '<':
             kargs.kss -= 0.1f;
             break;

         case 'L':
             kargs.kse += 0.1f;
             break;
         case '>':
             kargs.kse -= 0.1f;
             break;

        /*************************/
        /* METHODS               */
        /*************************/
        case '1':   //default
        case '2':   //classification
        case '3':   //bound
        case '4':   //depth color  
        case '5':   //tone
        case '6':   //mida
        case '7':   //sil
        case '8':   //tone+class
        case '9':   //dcb+class+light
            method = key - '0';
            break;
        
        
        case '`':
            fprintf(stdout, "got it\n");
            copyDepthBufferToHost(h_depthBuffer, d_depthBuffer, g_singleDimVolumeDataArrayLength);
            break;

        case '~':
            writeDBufToFile();
            break;

        default:
        	break;

    }

    printStuff();
    copyKParamsToDevice(&kargs, sizeof(kargs));
    glutPostRedisplay();
}

int ox, oy;
int buttonState = 0;

void mouse(int button, int state, int x, int y)
{
    if (state == GLUT_DOWN)
        buttonState  |= 1<<button;
    else if (state == GLUT_UP)
        buttonState = 0;

    ox = x; oy = y;
    glutPostRedisplay();
}

void motion(int x, int y)
{
    float dx, dy;
    dx = (float)(x - ox);
    dy = (float)(y - oy);

    if (buttonState == 4) {
        // right = zoom
        viewTranslation.z += dy / 100.0f;
    } 
    else if (buttonState == 2) {
        // middle = translate
        viewTranslation.x += dx / 100.0f;
        viewTranslation.y -= dy / 100.0f;
    }
    else if (buttonState == 1) {
        // left = rotate
        viewRotation.x += dy / 5.0f;
        viewRotation.y += dx / 5.0f;
    }

    ox = x; oy = y;
    glutPostRedisplay();
}

int iDivUp(int a, int b){
    return (a % b != 0) ? (a / b + 1) : (a / b);
}

void reshape(int w, int h)
{
    kargs.imageW = width = w; kargs.imageH = height = h;
    copyKParamsToDevice(&kargs, sizeof(kargs));
    
    initPixelBuffer();

    // calculate new grid size
    gridSize = dim3(iDivUp(width, blocksize.x), iDivUp(height, blocksize.y));



    glViewport(0, 0, w, h);

    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();

    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    glOrtho(0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
}

void cleanup()
{
    cutilCheckError( cutDeleteTimer( timer));

    freeCudaBuffers();

    if (pbo) {
        cudaGraphicsUnregisterResource(cuda_pbo_resource);
        glDeleteBuffersARB(1, &pbo);
        glDeleteTextures(1, &tex);
    }

    if (g_CheckRender) {
        delete g_CheckRender; g_CheckRender = NULL;
    }
}

void initGL(int *argc, char **argv)
{
    // initialize GLUT callback functions
    glutInit(argc, argv);
    glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE);
    glutInitWindowSize(width, height);
    glutCreateWindow("CUDA volume rendering");

    glewInit();
    if (!glewIsSupported("GL_VERSION_2_0 GL_ARB_pixel_buffer_object")) {
        shrLog("Required OpenGL extensions missing.");
	shrQAFinishExit(*argc, (const char **)argv, QA_WAIVED);
    }
}

void initPixelBuffer()
{
    if (pbo) {
        // unregister this buffer object from CUDA C
        cutilSafeCall(cudaGraphicsUnregisterResource(cuda_pbo_resource));

        // delete old buffer
        glDeleteBuffersARB(1, &pbo);
        glDeleteTextures(1, &tex);
    }

    // create pixel buffer object for display
    glGenBuffersARB(1, &pbo);
    glBindBufferARB(GL_PIXEL_UNPACK_BUFFER_ARB, pbo);
    glBufferDataARB(GL_PIXEL_UNPACK_BUFFER_ARB, width*height*sizeof(GLubyte)*4, 0, GL_STREAM_DRAW_ARB);
    glBindBufferARB(GL_PIXEL_UNPACK_BUFFER_ARB, 0);

    // register this buffer object with CUDA
	cutilSafeCall(cudaGraphicsGLRegisterBuffer(&cuda_pbo_resource, pbo, cudaGraphicsMapFlagsWriteDiscard));	

    // create texture for displayd_render
    glGenTextures(1, &tex);
    glBindTexture(GL_TEXTURE_2D, tex);
    glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA8, width, height, 0, GL_RGBA, GL_UNSIGNED_BYTE, NULL);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
    glBindTexture(GL_TEXTURE_2D, 0);
}

// Load raw data from disk
void *loadRawFile(char *filename, size_t size)
{
    FILE *fp = fopen(filename, "rb");
    if (!fp) {
        fprintf(stderr, "Error opening file '%s'\n", filename);
        return 0;
    }

    void *data = malloc(size);
    size_t read = fread(data, 1, size, fp);
    fclose(fp);

    shrLog("Read '%s', %d bytes\n", filename, read);

    return data;
}

// General initialization call for CUDA Device
int chooseCudaDevice(int argc, char **argv, bool bUseOpenGL)
{
    int result = 0;
    if (bUseOpenGL) {
        result = cutilChooseCudaGLDevice(argc, argv);
    } else {
        result = cutilChooseCudaDevice(argc, argv);
    }
    return result;
}

VolumeType 
getDataAtPt(VolumeType * volumedata, VolumeType x, VolumeType y, VolumeType z)
{
    x = Clamp(x, 0, volumeSize.width  - 1);
    y = Clamp(y, 0, volumeSize.height - 1);
    z = Clamp(z, 0, volumeSize.depth  - 1);
    
    return volumedata[x + (y * volumeSize.width) + (z * volumeSize.width * volumeSize.height)];
}

/*
 *	Calculate the gradient and gradient magnitude for each value.
 *  Uses central difference approximation, as explained by Levoy (1988).
 *
 *  volumedata: The volume data scalar values.
 *  gradMags:   The array to store the gradient magnitudes in.
 *  grads:      The array to store the gradients in. 
 *              Each gradient vector is stored as an ordered triple.
 */
void 
calculateGradient(VolumeType * volumedata, float * grads)
{
    int gradsIdx = 0;

    int depth  = (int)volumeSize.depth;
    int height = (int)volumeSize.height;
    int width  = (int)volumeSize.width;
    for (int z = 0; z < depth; z++)
    {
	    for (int y = 0; y < height; y++)
	    {
	        for (int x = 0; x < width; x++)
	        {
				VolumeType dxFor =  getDataAtPt( volumedata, x+1,y,z);
				VolumeType dxBack = getDataAtPt( volumedata, x-1,y,z);
				VolumeType dyFor =  getDataAtPt( volumedata, x,y+1,z);
				VolumeType dyBack = getDataAtPt( volumedata, x,y-1,z);
				VolumeType dzFor =  getDataAtPt( volumedata, x,y,z+1);
				VolumeType dzBack = getDataAtPt( volumedata, x,y,z-1);

				float dx = (dxBack - dxFor)/2.f;
				float dy = (dyBack - dyFor)/2.f;
				float dz = (dzBack - dzFor)/2.f;

				float mag = sqrtf(dx*dx + dy*dy + dz*dz);
				float len = (mag == 0) ? 1 : mag;
				grads[gradsIdx]   = dx/len;
				grads[gradsIdx+1] = dy/len;
				grads[gradsIdx+2] = dz/len;
				grads[gradsIdx+3] = mag;
				gradsIdx +=4;
	        }
	    }
    }
}

/***********************************
 * Initialize kargs to defaults
 ************************************/

void initKparams()
{
    method = g_options->find("m")->second;
    kargs.mida_gamma = 0.f;

	kargs.density = g_options->find("den")->second;
	kargs.brightness = g_options->find("brt")->second;
	kargs.transferOffset =  g_options->find("tranoff")->second;
	kargs.transferScale =  g_options->find("transcl")->second;
	kargs.imageW = width;
	kargs.imageH = height;

	kargs.lightx = g_options->find("ltx")->second;
    kargs.lighty = g_options->find("lty")->second;
    kargs.lightz = g_options->find("ltz")->second;

	kargs.fallOff = g_options->find("fo")->second;

	kargs.kgc = g_options->find("kgc")->second;
    kargs.kgs = g_options->find("kgs")->second; 
    kargs.kge = g_options->find("kge")->second;   //boundary

	kargs.ksc = g_options->find("ksc")->second;
    kargs.kss = g_options->find("kss")->second;
    kargs.kse = g_options->find("kse")->second;   //sil

	kargs.kds = g_options->find("kds")->second; 
    kargs.kde = g_options->find("kde")->second;                     //depth blend

	kargs.ktb = g_options->find("ktb")->second; 
    kargs.ktd = g_options->find("ktd")->second; 
    kargs.kty = g_options->find("kty")->second;   //tone

    kargs.mat_a = g_options->find("mat_a")->second; 
    kargs.mat_d = g_options->find("mat_d")->second; 
    kargs.mat_s = g_options->find("mat_s")->second;  //material values

    kargs.isoValues[0] = g_isoValues->at(0).isoval;
    kargs.isoValues[1] = g_isoValues->at(1).isoval;
    kargs.isoValues[2] = g_isoValues->at(2).isoval;
    kargs.isoAlphas[0] = g_isoValues->at(0).isoalpha;
    kargs.isoAlphas[1] = g_isoValues->at(1).isoalpha;
    kargs.isoAlphas[2] = g_isoValues->at(2).isoalpha;

    kargs.numIsoValues = 3;
}





////////////////////////////////////////////////////////////////////////////////
// Program main
////////////////////////////////////////////////////////////////////////////////
int
main( int argc, char** argv)
{
    pArgc = &argc;
    pArgv = argv;

    bool bTestResult = true;
    
    shrQAStart(argc, argv);

    //start logs
    shrSetLogFileName ("volumeRender.txt");
    shrLog("%s Starting...\n\n", argv[0]); 

    if (cutCheckCmdLineFlag(argc, (const char **)argv, "qatest") ||
        cutCheckCmdLineFlag(argc, (const char **)argv, "noprompt")) 
    {
        g_bQAReadback = true;
        fpsLimit = frameCheckNumber;
    }

    if (cutCheckCmdLineFlag(argc, (const char **)argv, "glverify")) 
    {
        g_bQAGLVerify = true;
        fpsLimit = frameCheckNumber;
    }

    if (g_bQAReadback) {
        // use command-line specified CUDA device, otherwise use device with highest Gflops/s
        chooseCudaDevice(argc, argv, false);
    } else {
        // First initialize OpenGL context, so we can properly set the GL for CUDA.
        // This is necessary in order to achieve optimal performance with OpenGL/CUDA interop.
        initGL( &argc, argv );

        // use command-line specified CUDA device, otherwise use device with highest Gflops/s
        chooseCudaDevice(argc, argv, true);
    }

    // parse arguments
    char *filename;
    if (cutGetCmdLineArgumentstr( argc, (const char**) argv, "file", &filename)) {
        volumeFilename = filename;
    }
    int n, n2;
    if (cutGetCmdLineArgumenti( argc, (const char**) argv, "size", &n)) {
        volumeSize.width = volumeSize.height = volumeSize.depth = n;
    }
    if (cutGetCmdLineArgumenti( argc, (const char**) argv, "xsize", &n)) {
        volumeSize.width = n;
    }
    if (cutGetCmdLineArgumenti( argc, (const char**) argv, "ysize", &n)) {
        volumeSize.height = n;
    }
    if (cutGetCmdLineArgumenti( argc, (const char**) argv, "zsize", &n)) {
         volumeSize.depth = n;
    }
    char *tfFile;
    if (cutGetCmdLineArgumentstr( argc, (const char**) argv, "trfile", &tfFile) &&
        cutGetCmdLineArgumenti  ( argc, (const char**) argv, "trlen",  &n) )  {
    	transFuncFileName = tfFile;
    	tfLines = n;
    }

//     if (cutGetCmdLineArgumenti(argc, (const char**)argv, "iv1", &n) &&
//         cutGetCmdLineArgumenti(argc, (const char**)argv, "ia1", &n2)) {
//             g_isoValues[0] = n;
//             g_isoAlphas[0] = n2;
//     }
//     if (cutGetCmdLineArgumenti(argc, (const char**)argv, "iv2", &n) &&
//         cutGetCmdLineArgumenti(argc, (const char**)argv, "ia2", &n2)) {
//             g_isoValues[1] = n;
//             g_isoAlphas[1] = n2;
//     }
//     if (cutGetCmdLineArgumenti(argc, (const char**)argv, "iv3", &n) &&
//         cutGetCmdLineArgumenti(argc, (const char**)argv, "ia3", &n2)) {
//             g_isoValues[2] = n;
//             g_isoAlphas[2] = n2;
//     }
//     if (cutGetCmdLineArgumenti( argc, (const char**) argv, "m", &n)) {
//         g_method = n;
//     }
//     float pos;
//     if (cutGetCmdLineArgumentf( argc, (const char**) argv, "lightx", &pos)) {
//         g_lightpos.x = pos;
//     }
//     if (cutGetCmdLineArgumentf( argc, (const char**) argv, "lighty", &pos)) {
//         g_lightpos.y = pos;
//     }
//     if (cutGetCmdLineArgumentf( argc, (const char**) argv, "lightz", &pos)) {
//         g_lightpos.z = pos;
//     }
//     if (cutGetCmdLineArgumentf( argc, (const char**) argv, "mg", &pos)) {
//         g_mida_gama = pos;
//     }

  

   
 
    // load volume data
    char* path = shrFindFilePath(volumeFilename, argv[0]);
    if (path == 0) {
        shrLog("Error finding file '%s'\n", volumeFilename);
        shrQAFinishExit(argc, (const char **)argv, QA_FAILED);
    }
    
    //load transfer function
	float4 *h_transFunc =  (float4*)malloc(sizeof(float4)*tfLines);
    g_isoValues = new std::vector<IsoPair>();
    g_transValues = new std::vector<float4>();
    g_options = new std::map<string,float>();
   
	int tfLinesParsed = parseTrans(transFuncFileName, g_options, g_transValues, g_isoValues);
	if (tfLinesParsed != tfLines){
		shrLog("Error importing transfer function: Incorrect number of lines given (%d) vs. lines found (%d)\n Exiting...\n", tfLines, tfLinesParsed);
		return 0;
	}
    
    for (int i = 0; i < g_transValues->size(); i++)
    {
        h_transFunc[i] = g_transValues->at(i);
    }



    size_t length = volumeSize.width*volumeSize.height*volumeSize.depth;
    g_singleDimVolumeDataArrayLength = length;
    size_t size = length*sizeof(VolumeType);
    void *h_volume = loadRawFile(path, size);
    float *h_grads = (float*) malloc(sizeof(float) * length*4); //ordered quads
    

    calculateGradient((VolumeType*)h_volume, h_grads);
    initCuda(h_volume, h_grads, h_transFunc, NULL, tfLines,volumeSize);
    cutilSafeCall( cudaMalloc((void**)&d_depthBuffer, volumeSize.width*volumeSize.height*volumeSize.depth) );
    //cutilSafeCall( cudaMemset((void**)&d_depthBuffer, 0, volumeSize.width*volumeSize.height) );
    
    initKparams();
    
    free(h_volume);
    free(h_grads);
    free(h_transFunc);

    h_depthBuffer = (float *)malloc(sizeof(float) * length);

    cutilCheckError( cutCreateTimer( &timer ) );
	cutilCheckError(cutCreateTimer(&fpsTimer));

    shrLog("Press '+' and '-' to change density (0.01 increments)\n"
           "      '9' and '0' to change brightness\n"
           "      '[' and ']' to modify transfer function offset\n"
           "      'o' and 'p' to modify transfer function scale\n\n");

    // calculate new grid size
    gridSize = dim3(iDivUp(width, blocksize.x), iDivUp(height, blocksize.y));
    
    printStuff();

//     //if (g_bQAReadback) {
//  	if (false){
//          g_CheckRender = new CheckBackBuffer(width, height, 4, false);
//          g_CheckRender->setPixelFormat(GL_RGBA);
//          g_CheckRender->setExecPath(argv[0]);
//          g_CheckRender->EnableQAReadback(true);
//  
//          uint *d_output;
//          cutilSafeCall(cudaMalloc((void**)&(d_output), width*height*sizeof(uint)));
//          cutilSafeCall(cudaMemset(d_output, 0, width*height*sizeof(uint)));
//  
//          float modelView[16] = 
//          {
//              1.0f, 0.0f, 0.0f, 0.0f,
//              0.0f, 1.0f, 0.0f, 0.0f,
//              0.0f, 0.0f, 1.0f, 0.0f,
//              0.0f, 0.0f, 4.0f, 1.0f
//          };
//  
//          invViewMatrix[0] = modelView[0]; invViewMatrix[1] = modelView[4]; invViewMatrix[2] = modelView[8]; invViewMatrix[3] = modelView[12];
//          invViewMatrix[4] = modelView[1]; invViewMatrix[5] = modelView[5]; invViewMatrix[6] = modelView[9]; invViewMatrix[7] = modelView[13];
//          invViewMatrix[8] = modelView[2]; invViewMatrix[9] = modelView[6]; invViewMatrix[10] = modelView[10]; invViewMatrix[11] = modelView[14];
//  
//          // call CUDA kernel, writing results to PBO
//  	    copyInvViewMatrix(invViewMatrix, sizeof(float4)*3);
//          
//          // Start timer 0 and process n loops on the GPU 
//          int nIter = 10;
//          for (int i = -1; i < nIter; i++)
//          {
//              if( i == 0 ) {
//                  cutilDeviceSynchronize();
//                  cutStartTimer(timer); 
//              }
//              
//              render_kernel(gridSize, blocksize, d_output, d_depthBuffer, method);
//          }
//          cutilDeviceSynchronize();
//          cutStopTimer(timer);
//          // Get elapsed time and throughput, then log to sample and master logs
//          double dAvgTime = cutGetTimerValue(timer)/(nIter * 1000.0);
//          shrLogEx(LOGBOTH | MASTER, 0, "volumeRender, Throughput = %.4f MTexels/s, Time = %.5f s, Size = %u Texels, NumDevsUsed = %u, Workgroup = %u\n", 
//                 (1.0e-6 * width * height)/dAvgTime, dAvgTime, (width * height), 1, blocksize.x * blocksize.y);
//          
//  
//          cutilCheckMsg("Error: render_kernel() execution FAILED");
//          cutilSafeCall( cutilDeviceSynchronize() );
//  
//          cutilSafeCall( cudaMemcpy(g_CheckRender->imageData(), d_output, width*height*4, cudaMemcpyDeviceToHost) );
//          g_CheckRender->savePPM(sOriginal[g_Index], true, NULL);
//  
//          bTestResult = g_CheckRender->PPMvsPPM(sOriginal[g_Index], sReference[g_Index], MAX_EPSILON_ERROR, THRESHOLD);
//  
//          cudaFree(d_output);
//      	freeCudaBuffers();
//  
//          if (g_CheckRender) {
//              delete g_CheckRender; g_CheckRender = NULL;
//          }
//  
//      } else {
        // This is the normal rendering path for VolumeRender
        glutDisplayFunc(display);
        glutKeyboardFunc(keyboard);
        glutMouseFunc(mouse);
        glutMotionFunc(motion);
        glutReshapeFunc(reshape);
        glutIdleFunc(idle);

        initPixelBuffer();

        if (g_bQAGLVerify) {
            g_CheckRender = new CheckBackBuffer(width, height, 4);
            g_CheckRender->setPixelFormat(GL_RGBA);
            g_CheckRender->setExecPath(argv[0]);
            g_CheckRender->EnableQAReadback(true);
        }
        atexit(cleanup);

        copyKParamsToDevice(&kargs, sizeof(kargs));
		cutStartTimer(fpsTimer);
		glutMainLoop();
    //}
    cutilDeviceReset();
    shrQAFinishExit(argc, (const char **)argv, (bTestResult ? QA_PASSED : QA_FAILED) );

   
}

//TODO: double check col_l uses the right equation (missing anything?).
//TODO: check grads against vtk, or something that calcs gradients.
//TODO: implement transfer function technique from Levoy 1988.
//TODO: convolute the depth buffer?!?!
//TODO: comment method functions with sources of each method.
//TODO: improve gradient calculations.
//TODO: free d_depthBuffer, h_depthBuffer
